package airportpainter.earth;

/*  class Magfield.java */

/* Module to calculate magnetic variation and field given position,
**               altitude, and date
** Implements the NIMA (formerly DMA) WMM and IGRF models
**
**    http://www.nima.mil/GandG/ngdc-wmm2000.html
**    For WMM2000 coefficients:
**    ftp://ftp.ngdc.noaa.gov/Solid_Earth/Mainfld_Mag/DoD_Model/wmm.cof
**    For IGRF/DGRF coefficients:
**    http://swdcdb.kugi.kyoto-u.ac.jp/igrf/coef/igrfall.d
**
** Copyright (C) 2000  Edward A Williams <Ed_Williams@compuserve.com>
**
**  The routine uses a spherical harmonic expansion of the magnetic
** potential up to twelfth order, together with its time variation, as
** described in Chapter 4 of "Geomagnetism, Vol 1, Ed. J.A.Jacobs,
** Academic Press (London 1987)". The program first converts geodetic
** coordinates (lat/long on elliptic earth and altitude) to spherical
** geocentric (spherical lat/long and radius) coordinates. Using this,
** the spherical (B_r, B_theta, B_phi) magnetic field components are
** computed from the model. These are finally referred to surface (X, Y,
** Z) coordinates.
**
**   Fields are accurate to better than 200nT, variation and dip to
** better than 0.5 degrees, with the exception of the declination near
** the magnetic poles (where it is ill-defined) where the error may reach
** 4 degrees or more.
**
**   Variation is undefined at both the geographic and
** magnetic poles, even though the field itself is well-behaved. To
** avoid the routine blowing up, latitude entries corresponding to
** the geographic poles are slightly offset. At the magnetic poles,
** the routine returns zero variation.
**
** HISTORY
** Adapted from EAW Excel 3.0 version 3/27/94 EAW
** Recoded in C++ by Starry Chan
** WMM95 added and rearranged in ANSI-C EAW 7/9/95
** Put shell around program and made Borland & GCC compatible EAW 11/22/95
** IGRF95 added 2/96 EAW
** WMM2000 IGR2000 added 2/00 EAW
** Released under GPL  3/26/00 EAW
** Adaptions and modifications for the SimGear project  3/27/2000 CLO
** Removed all pow() calls and made static roots[][] arrays to
** save many sqrt() calls on subsequent invocations
** 3/28/2000  Norman Vine -- nhv@yahoo.com
** Put in some bullet-proofing to handle magnetic and geographic poles.
** 3/28/2000 EAW
** Converted to Java class
** 12/6/2000 Reece Robinson
**
** This program is free software; you can redistribute it and/or
** modify it under the terms of the GNU General Public License as
** published by the Free Software Foundation; either version 2 of the
** License, or (at your option) any later version.
**
** This program is distributed in the hope that it will be useful, but
** WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
** General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**
*/
import java.lang.Math;
import java.util.Calendar;
import java.util.GregorianCalendar;

public class Magfield {

  private static double pi = 3.14159265358979;
  private static final double a = 6378.16;	/* major radius (km) IAU66 ellipsoid */
//  private static final double f = 1.0 / 298.25;	/* inverse flattening IAU66 ellipsoid */
  private static final double b = 6378.16 * (1.0 -1.0 / 298.25 );
	/* minor radius b=a*(1-f) */
  private static final double r_0 = 6371.2;	/* "mean radius" for spherical harmonic expansion */

/*IGRF90 constants */
  private static final double[][] gnm_igrf90 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29775.4, -1851, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-2135.8, 3058.2, 1693.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{1314.6, -2240.2, 1245.6, 806.5, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{938.9, 782.3, 323.9, -422.7, 141.7, 0, 0, 0, 0, 0, 0, 0, 0},
	{-211, 352.5, 243.8, -110.8, -165.6, -37, 0, 0, 0, 0, 0, 0, 0},
	{60.7, 63.9, 60.4, -177.5, 2, 16.7, -96.3, 0, 0, 0, 0, 0, 0},
	{76.6, -64.2, 3.7, 27.5, 0.9, 5.7, 9.8, -0.5, 0, 0, 0, 0, 0},
	{22.4, 5.1, -0.9, -10.8, -12.4, 3.8, 3.8, 2.6, -6, 0, 0, 0, 0},
	{4.4, 9.9, 0.8, -12, 9.3, -3.9, -1.4, 7.3, 1.5, -5.5, 0, 0, 0},
	{-3.6, -3.9, 2.4, -5.3, -2.4, 4.4, 3, 1.2, 2.2, 2.9, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
  private static final double[][] hnm_igrf90 =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 5410.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -2277.7, -380, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -286.5, 293.3, -348.5, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 248.1, -239.5, 87, -299.4, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 47.2, 153.5, -154.4, -69.2, 97.7, 0, 0, 0, 0, 0, 0, 0},
	{0, -15.8, 82.7, 68.3, -52.5, 1.8, 26.9, 0, 0, 0, 0, 0, 0},
	{0, -81.1, -27.3, 0.6, 20.4, 16.4, -22.6, -5, 0, 0, 0, 0, 0},
	{0, 9.7, -19.9, 7.1, -22.1, 11.9, 11, -16, -10.7, 0, 0, 0, 0},
	{0, -20.8, 15.4, 9.5, -5.7, -6.4, 8.6, 9.1, -6.6, 1.9, 0, 0, 0},
	{0, 1.3, 0.4, 3.1, 5.6, -4.2, -0.5, -1.5, 3.8, -0.5, -6.2, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
  private static final double[][] gtnm_igrf90 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{18, 10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-12.9, 2.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{3.3, -6.7, 0, -5.9, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0.5, 0.6, -7, 0.5, -5.5, 0, 0, 0, 0, 0, 0, 0, 0},
	{0.6, -0.1, -1.6, -3.1, 0, 2.3, 0, 0, 0, 0, 0, 0, 0},
	{1.3, -0.2, 1.8, 1.3, -0.2, 0.1, 1.2, 0, 0, 0, 0, 0, 0},
	{0.6, -0.5, -0.3, 0.6, 1.6, 0.2, 0.2, 0.3, 0, 0, 0, 0, 0},
	{0.2, -0.7, -0.2, 0.1, -1.1, 0, 0, -0.5, -0.6, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
  private static final double[][] htnm_igrf90 =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -16.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -15.8, -13.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 4.4, 1.6, -10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 2.6, 1.8, 3.1, -1.4, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -0.1, 0.5, 0.4, 1.7, 0.4, 0, 0, 0, 0, 0, 0, 0},
	{0, 0.2, -1.3, 0, -0.9, 0.5, 1.2, 0, 0, 0, 0, 0, 0},
	{0, 0.6, 0.2, 0.8, -0.5, -0.2, 0, 0, 0, 0, 0, 0, 0},
	{0, 0.5, -0.2, 0.3, 0.3, 0.4, -0.5, -0.3, 0.6, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
/* IGRF95 */
  private static final double[][] gnm_igrf95 =
{
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {-29682.0, -1789.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {-2197.0, 3074.0, 1685.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {1329.0, -2268.0, 1249.0, 769.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {941.0, 782.0, 291.0, -421.0, 116.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {-210.0, 352.0, 237.0, -122.0, -167.0, -26.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0},
  {66.0, 64.0, 65.0, -172.0, 2.0, 17.0, -94.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {78.0, -67.0, 1.0, 29.0, 4.0, 8.0, 10.0, -2.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {24.0, 4.0, -1.0, -9.0, -14.0, 4.0, 5.0, 0.0, -7.0, 0.0,0.0, 0.0, 0.0},
  {4.0, 9.0, 1.0, -12.0, 9.0, -4.0, -2.0, 7.0, 0.0, -6.0, 0.0, 0.0, 0.0},
  {-3.0, -4.0, 2.0, -5.0, -2.0, 4.0, 3.0, 1.0, 3.0, 3.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0}
 };

  private static final double[][] hnm_igrf95 =
{
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 5318.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -2356.0, -425.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -263.0, 302.0, -406.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 262.0, -232.0, 98.0, -301.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 44.0, 157.0, -152.0, -64.0, 99.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -16.0, 77.0, 67.0, -57.0, 4.0, 28.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -77.0, -25.0, 3.0, 22.0, 16.0, -23.0, -3.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 12.0, -20.0, 7.0, -21.0, 12.0, 10.0, -17.0, -10.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -19.0, 15.0, 11.0, -7.0, -7.0, 9.0, 7.0, -8.0, 1.0, 0.0, 0.0, 0.0},
  {0.0, 2.0, 1.0, 3.0, 6.0, -4.0, 0.0, -2.0, 3.0, -1.0, -6.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0}
};

  private static final double[][] gtnm_igrf95 =
{
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {17.6, 13.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {-13.2, 3.7, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {1.5, -6.4, -0.2, -8.1, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.8, 0.9, -6.9, 0.5, -4.6, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.8, 0.1, -1.5, -2.0, -0.1, 2.3, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.5, -0.4, 0.6, 1.9, -0.2, -0.2, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {-0.2, -0.8, -0.6, 0.6, 1.2, 0.1, 0.2, -0.6, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.3, -0.2, 0.1, 0.4, -1.1, 0.3, 0.2, -0.9, -0.3, 0.0,0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0}
 };

  private static final double[][] htnm_igrf95 =
{
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -18.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -15.0, -8.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 4.1, 2.2, -12.1, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 1.8, 1.2, 2.7, -1.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.2, 1.2, 0.3, 1.8, 0.9, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.3, -1.6, -0.2, -0.9, 1.0, 2.2, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.8, 0.2, 0.6, -0.4, 0.0, -0.3, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.4, -0.2, 0.2, 0.7, 0.0, -1.2, -0.7, -0.6, 0.0,0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0}
};

/*WMM85 constants */
  private static final double[][] gnm_wmm85 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29879.8, -1903.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-2070.6, 3040.8, 1696.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{1303.9, -2203, 1241.7, 839.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{933.8, 781.8, 359, -424.5, 164.5, 0, 0, 0, 0, 0, 0, 0, 0},
	{-216.4, 353, 254.3, -93.7, -157.5, -45.2, 0, 0, 0, 0, 0, 0, 0},
	{53.2, 63.8, 51.3, -188.4, 3.3, 20.3, -101.7, 0, 0, 0, 0, 0, 0},
	{76.9, -60.7, 0.7, 25.4, -8.1, 6.9, 7, -4.4, 0, 0, 0, 0, 0},
	{18.4, 5.1, 1.2, -12, -9.1, 0.1, 4.7, 6.5, -9.5, 0, 0, 0, 0},
	{5.7, 10.9, 0.9, -12.2, 9.5, -3.3, -1, 6.5, 1.5, -4.8, 0, 0, 0},
	{-3.4, -4.7, 2.5, -5.5, -2.1, 4.6, 3.2, 0.6, 1.9, 2.8, -0.2, 0, 0},
	{2.3, -0.8, -2, 2.1, 0.2, -0.4, -0.4, 1.6, 1.5, -0.7, 2.3, 3.5, 0},
	{-1.8, 0, 0.1, -0.3, 0.5, 0.5, -0.6, -0.4, 0, -0.5, 0, 0.7, -0.2}};
  private static final double[][] hnm_wmm85 =
	{
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 5490.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -2189.1, -312, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -310.3, 282.6, -299.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 227.2, -246.7, 72.5, -299.1, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 43.4, 148.2, -154.8, -71.8, 91.5, 0, 0, 0, 0, 0, 0, 0},
	{0, -12.3, 87.9, 67.8, -51.1, -4, 20.8, 0, 0, 0, 0, 0, 0},
	{0, -80.1, -25.9, -0.9, 21.6, 18.5, -20, -7.7, 0, 0, 0, 0, 0},
	{0, 3.8, -20.2, 5, -24.2, 12.2, 7.6, -16.3, -10.9, 0, 0, 0, 0},
	{0, -20.8, 15.8, 9, -5, -6.4, 9.1, 9.9, -5.8, 2.3, 0, 0, 0},
	{0, 1.2, 0.4, 2.5, 5.6, -4.4, -0.5, -1.6, 3.7, -0.5, -6.1, 0, 0},
	{0, 1.3, 2, -1.1, -2.8, 0.7, -0.1, -2.4, -0.4, -1.5, -1.5, 0.7, 0},
	{0, 0.3, 0.6, 2.5, -1.7, 0.3, 0.2, -0.1, 0.1, 0.1, -1.4, 0.4, 0.7}};
  private static final double[][] gtnm_wmm85 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{21.9, 10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-11.2, 1.8, 9.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{8.3, -2, -0.6, 2.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-1.2, 0.1, -9.7, -1.7, -9.3, 0, 0, 0, 0, 0, 0, 0, 0},
	{1.4, -0.5, -1.2, -2.2, 0.9, 0, 0, 0, 0, 0, 0, 0, 0},
	{3.1, 0, 1.8, -0.2, -0.4, 2.4, 1.8, 0, 0, 0, 0, 0, 0},
	{-0.1, -0.8, -1.2, 1.1, 0, 0.6, -1.8, -1.2, 0, 0, 0, 0, 0},
	{0.2, 0, 0.7, 0.1, 0.2, -0.3, -0.1, 0.2, -2.2, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
  private static final double[][] htnm_wmm85 =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -31.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -9.7, -19.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 6.1, 1.3, -13, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 1.3, 3.6, 2.5, 0.6, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -0.9, 0.6, 0.3, 2.4, -1.4, 0, 0, 0, 0, 0, 0, 0},
	{0, 0.7, -2.1, -1.4, -4.3, -0.7, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 1.2, 2, 2.6, 0.9, 0.8, 0.4, 0, 0, 0, 0, 0},
	{0, -0.6, -1.5, 0.1, -1.1, 0.4, -2, 0.9, 1.5, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};

/* wmm90 constants */
  private static final double[][] gnm_wmm90 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29780.5, -1851.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-2134.3, 3062.2, 1691.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{1312.9, -2244.7, 1246.8, 808.6, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{933.5, 784.9, 323.5, -421.7, 139.2, 0, 0, 0, 0, 0, 0, 0, 0},
	{-208.3, 352.2, 246.5, -110.8, -162.3, -37.2, 0, 0, 0, 0, 0, 0, 0},
	{59, 63.7, 60, -181.3, 0.4, 15.4, -96, 0, 0, 0, 0, 0, 0},
	{76.1, -62.1, 1.3, 30.2, 4.7, 7.9, 10.1, 1.9, 0, 0, 0, 0, 0},
	{22.9, 2.3, -1.2, -11.7, -17.5, 2.2, 5.7, 3, -7, 0, 0, 0, 0},
	{3.6, 9.5, -0.9, -10.7, 10.7, -3.2, -1.4, 6.3, 0.8, -5.5, 0, 0, 0},
	{-3.3, -2.6, 4.5, -5.6, -3.6, 3.9, 3.2, 1.7, 3, 3.7, 0.7, 0, 0},
	{1.3, -1.4, -2.5, 3.2, 0.2, -1.1, 0.3, -0.3, 0.9, -1.1, 2.4, 3, 0},
	{-1.3, 0.1, 0.5, 0.7, 0.4, -0.2, -1.1, 0.9, -0.6, 0.8, 0.2, 0.4, 0.2}};
  private static final double[][] hnm_wmm90 =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 5407.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -2278.3, -384.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -284.9, 291.7, -352.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 249.4, -232.7, 91.3, -296.5, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 40.8, 148.7, -154.6, -67.6, 97.4, 0, 0, 0, 0, 0, 0, 0},
	{0, -14.7, 82.2, 70, -56.2, -1.4, 24.6, 0, 0, 0, 0, 0, 0},
	{0, -78.6, -26.7, 0.1, 19.9, 17.9, -21.5, -6.8, 0, 0, 0, 0, 0},
	{0, 9.7, -19.3, 6.6, -20.1, 13.4, 9.8, -19, -9.1, 0, 0, 0, 0},
	{0, -21.9, 14.3, 9.5, -6.7, -6.4, 9.1, 8.9, -8, 2.1, 0, 0, 0},
	{0, 2.6, 1.2, 2.6, 5.7, -4, -0.4, -1.7, 3.8, -0.8, -6.5, 0, 0},
	{0, 0, 1, -1.6, -2.2, 1.1, -0.7, -1.7, -1.5, -1.3, -1.1, 0.6, 0},
	{0, 0.7, 0.7, 1.3, -1.5, 0.3, 0.2, -1.1, 1.2, -0.2, -1.3, 0.6, 0.6}};
  private static final double[][] gtnm_wmm90 =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{16, 9.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-11.7, 3.7, 1.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{2.1, -7.6, 0, -5.8, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-0.8, 1, -7.4, 0.8, -6.4, 0, 0, 0, 0, 0, 0, 0, 0},
	{1.7, 0, 0, -2.7, 0, 3, 0, 0, 0, 0, 0, 0, 0},
	{0.8, 0, 1.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0.5, 0, -0.9, 1.5, 2.7, -1, 0, 0, 0, 0, 0, 0, 0},
	{0, -1.1, 0, 0, -2.1, 0, 1, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
  private static final double[][] htnm_wmm90 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -13.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -12.8, -14.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 3.1, 0.8, -11.3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 3.3, 3.7, 2.8, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, -2.1, 1.2, 1.2, 0.6, 0, 0, 0, 0, 0, 0, 0},
	{0, -0.6, -0.6, 0, -2.3, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0.6, 0.8, 0, 0, 0, 0.4, 0, 0, 0, 0, 0, 0},
	{0, 0.4, -0.8, 0.5, 0.3, 0.5, 0, -0.7, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
/* wmm95 constants */
  private static final double[][] gnm_wmm95 =
  {
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  {-29682.1,-1782.2,0,0,0,0,0,0,0,0,0,0,0},
  {-2194.7,3078.6,1685.7,0,0,0,0,0,0,0,0,0,0},
  {1318.8,-2273.6,1246.9,766.3,0,0,0,0,0,0,0,0,0},
  {940,782.9,290.9,-418.9,113.8,0,0,0,0,0,0,0,0},
  {-209.5,354,238.2,-122.1,-162.8,-23.3,0,0,0,0,0,0,0},
  {68.5,65.6,64.1,-169.1,-0.5,16.5,-91,0,0,0,0,0,0},
  {78,-68.1,0.1,29.6,6,8.7,9.2,-2.4,0,0,0,0,0},
  {24.7,3.4,-1.5,-9.6,-16.5,2.6,3.6,-4.9,-8.5,0,0,0,0},
  {2.9,7.5,0.4,-10.3,9.7,-2.3,-2.4,6.8,-0.5,-6.5,0,0,0},
  {-2.9,-3.3,2.8,-4.3,-3.1,2.4,2.8,0.7,4.1,3.6,0.6,0,0},
  {1.7,-1.6,-3.6,1.2,-0.6,0.1,-0.7,-0.8,1.3,-0.3,2.2,4.2,0},
  {-1.8,0.9,-0.1,-0.5,0.8,0.2,0.5,0.4,-0.4,0.3,0.2,0.4,0.6}};
  private static final double[][] hnm_wmm95 =
	{
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  {0,5315.6,0,0,0,0,0,0,0,0,0,0,0},
  {0,-2359.1,-418.6,0,0,0,0,0,0,0,0,0,0},
  {0,-261.1,301,-416.5,0,0,0,0,0,0,0,0,0},
  {0,259.4,-230.9,99.8,-306.1,0,0,0,0,0,0,0,0},
  {0,43.7,157.6,-150.1,-59.2,104.4,0,0,0,0,0,0,0},
  {0,-15.2,74.3,69.4,-55.3,3,33.3,0,0,0,0,0,0},
  {0,-76.1,-24.5,1.6,20,16.5,-23.6,-6.8,0,0,0,0,0},
  {0,14.9,-19.5,6.3,-20.4,12.2,7,-19,-8.8,0,0,0,0},
  {0,-19.8,14.6,10.9,-7.5,-6.8,9.3,7.7,-8.1,2.6,0,0,0},
  {0,3.2,1.7,2.9,5.6,-3.4,-0.7,-2.9,2.3,-1.6,-6.6,0,0},
  {0,0.3,1,-3.6,-1.4,1.9,0.2,-1.3,-2.4,-0.6,-2.2,1.3,0},
  {0,0.3,1.4,0.8,-3,0.7,0.5,-0.8,0.6,0.1,-1.3,-0.4,0.9}};
  private static final double[][] gtnm_wmm95 =
  {
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  {17.6,13.2,0,0,0,0,0,0,0,0,0,0,0},
  {-13.7,4,-0.3,0,0,0,0,0,0,0,0,0,0},
  {0.8,-6.6,-0.5,-8.5,0,0,0,0,0,0,0,0,0},
  {1.2,1.1,-6.8,0.3,-4.5,0,0,0,0,0,0,0,0},
  {0.9,0.5,-1.4,-1.7,0,2.1,0,0,0,0,0,0,0},
  {0.4,-0.3,0.3,2.1,0,-0.4,-0.4,0,0,0,0,0,0},
  {-0.3,-1.1,-0.5,0.5,1.3,0.1,0,-0.9,0,0,0,0,0},
  {0.1,0,0.4,0.3,-1.3,0.5,0.4,-0.9,0.1,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0}};
  private static final double[][] htnm_wmm95 =
  {
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  {0,-18,0,0,0,0,0,0,0,0,0,0,0},
  {0,-14.6,-7.2,0,0,0,0,0,0,0,0,0,0},
  {0,4,2.2,-12.6,0,0,0,0,0,0,0,0,0},
  {0,1.3,1,2.5,-1.2,0,0,0,0,0,0,0,0},
  {0,0.5,1.5,0.6,1.7,0.6,0,0,0,0,0,0,0},
  {0,0.7,-1.5,-0.5,-0.7,1.1,2.6,0,0,0,0,0,0},
  {0,0.3,0,0.7,-0.6,0.1,-0.6,-0.4,0,0,0,0,0},
  {0,0.4,-0.3,0.1,0.8,-0.1,-1.3,-0.9,-1.1,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0},
  {0,0,0,0,0,0,0,0,0,0,0,0,0}};

  private static final double[][] gnm_wmm2000 =
        {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-29616.0, -1722.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-2266.7, 3070.2, 1677.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1322.4, -2291.5, 1255.9, 724.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{932.1, 786.3, 250.6, -401.5, 106.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-211.9, 351.6, 220.8, -134.5, -168.8, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{73.8, 68.2, 74.1, -163.5, -3.8, 17.1, -85.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{77.4, -73.9, 2.2, 35.7, 7.3, 5.2, 8.4, -1.5, 0.0, 0.0, 0.0, 0.0, 0.0},
	{23.3, 7.3, -8.5, -6.6, -16.9, 8.6, 4.9, -7.8, -7.6, 0.0, 0.0, 0.0, 0.0},
	{5.7, 8.5, 2.0, -9.8, 7.6, -7.0, -2.0, 9.2, -2.2, -6.6, 0.0, 0.0, 0.0},
	{-2.2, -5.7, 1.6, -3.7, -0.6, 4.1, 2.2, 2.2, 4.6, 2.3, 0.1, 0.0, 0.0},
	{3.3, -1.1, -2.4, 2.6, -1.3, -1.7, -0.6, 0.4, 0.7, -0.3, 2.3, 4.2, 0.0},
	{-1.5, -0.2, -0.3, 0.5, 0.2, 0.9, -1.4, 0.6, -0.6, -1.0, -0.3, 0.3, 0.4},
	};

  private static final double[][] hnm_wmm2000 =
        {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 5194.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -2484.8, -467.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -224.7, 293.0, -486.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 273.3, -227.9, 120.9, -302.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 42.0, 173.8, -135.0, -38.6, 105.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -17.4, 61.2, 63.2, -62.9, 0.2, 43.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -62.3, -24.5, 8.9, 23.4, 15.0, -27.6, -7.8, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 12.4, -20.8, 8.4, -21.2, 15.5, 9.1, -15.5, -5.4, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.4, 13.9, 12.0, -6.2, -8.6, 9.4, 5.0, -8.4, 3.2, 0.0, 0.0, 0.0},
	{0.0, 0.9, -0.7, 3.9, 4.8, -5.3, -1.0, -2.4, 1.3, -2.3, -6.4, 0.0, 0.0},
	{0.0, -1.5, 0.7, -1.1, -2.3, 1.3, -0.6, -2.8, -1.6, -0.1, -1.9, 1.4, 0.0},
	{0.0, -1.0, 0.7, 2.2, -2.5, -0.2, 0.0, -0.2, 0.0, 0.2, -0.9, -0.2, 1.0},
	};

  private static final double[][] gtnm_wmm2000 =
        {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{14.7, 11.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-13.6, -0.7, -1.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.3, -4.3, 0.9, -8.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-1.6, 0.9, -7.6, 2.2, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.9, -0.2, -2.5, -2.7, -0.9, 1.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1.2, 0.2, 1.7, 1.6, -0.1, -0.3, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.4, -0.8, -0.2, 1.1, 0.4, 0.0, -0.2, -0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.3, 0.6, -0.8, 0.3, -0.2, 0.5, 0.0, -0.6, 0.1, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

  private static final double[][] htnm_wmm2000 =
  {
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -21.5, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 6.4, -1.3, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 2.3, 0.7, 3.7, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 2.1, 2.3, 3.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.3, -1.7, -0.9, -1.0, -0.1, 1.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 1.4, 0.2, 0.7, 0.4, -0.3, -0.8, -0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.5, 0.1, -0.2, 0.0, 0.1, -0.1, 0.3, 0.2, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  };

  private static final double[][] gnm_igrf2000 =
  {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-29615.0, -1728.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-2267.0, 3072.0, 1672.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1341.0, -2290.0, 1253.0, 715.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{935.0, 787.0, 251.0, -405.0, 110.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-217.0, 351.0, 222.0, -131.0, -169.0, -12.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{72.0, 68.0, 74.0, -161.0, -5.0, 17.0, -91.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{79.0, -74.0, 0.0, 33.0, 9.0, 7.0, 8.0, -2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{25.0, 6.0, -9.0, -8.0, -17.0, 9.0, 7.0, -8.0, -7.0, 0.0, 0.0, 0.0, 0.0},
	{5.0, 9.0, 3.0, -8.0, 6.0, -9.0, -2.0, 9.0, -4.0, -8.0, 0.0, 0.0, 0.0},
	{-2.0, -6.0, 2.0, -3.0, 0.0, 4.0, 1.0, 2.0, 4.0, 0.0, -1.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

  private static final double[][] hnm_igrf2000 =
  {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 5186.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -2478.0, -458.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -227.0, 296.0, -492.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 272.0, -232.0, 119.0, -304.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 44.0, 172.0, -134.0, -40.0, 107.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -17.0, 64.0, 65.0, -61.0, 1.0, 44.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -65.0, -24.0, 6.0, 24.0, 15.0, -25.0, -6.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 12.0, -22.0, 8.0, -21.0, 15.0, 9.0, -16.0, -3.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.0, 13.0, 12.0, -6.0, -8.0, 9.0, 4.0, -8.0, 5.0, 0.0, 0.0, 0.0},
	{0.0, 1.0, 0.0, 4.0, 5.0, -6.0, -1.0, -3.0, 0.0, -2.0, -8.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

  private static final double[][] gtnm_igrf2000 =
  {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{14.6, 10.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-12.4, 1.1, -1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.7, -5.4, 0.9, -7.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-1.3, 1.6, -7.3, 2.9, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.7, -2.1, -2.8, -0.8, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1.0, -0.4, 0.9, 2.0, -0.6, -0.3, 1.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.4, -0.4, -0.3, 1.1, 1.1, -0.2, 0.6, -0.9, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.3, 0.2, -0.3, 0.4, -1.0, 0.3, -0.5, -0.7, -0.4, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

  private static final double[][] htnm_igrf2000 =
  {
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -22.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.6, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 6.0, -0.1, -14.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 2.1, 1.3, 5.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.1, 0.6, 1.7, 1.9, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.2, -1.4, 0.0, -0.8, 0.0, 0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 1.1, 0.0, 0.3, -0.1, -0.6, -0.7, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.1, 0.0, 0.0, 0.3, 0.6, -0.4, 0.3, 0.7, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

  private static final int nmax = 12;

  private static double[][] P = new double[13][13];
  private static double[][] DP = new double[13][13];
  private static double[][] gnm = new double[13][13];
  private static double[][] hnm = new double[13][13];
  private static double[] sm = new double[13];
  private static double[] cm = new double[13];

  private static double[] root = new double[13];
  private static double[][][] roots = new double[13][13][2];

/* Convert date to Julian day */
  public static long yymmdd_to_julian_days(int yy,int mm,int dd) {
    long jd;

    yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
    jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
    jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
    jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;

    return(jd);
  }

/* Convert degrees to radians */
  public static double deg_to_rad(double deg) {
    return Math.toRadians(deg);
  }

/* Convert radians to degrees */
  public static double rad_to_deg(double rad) {
    return Math.toDegrees(rad);
  }

/*
 * return variation (in radians) given geodetic latitude (radians), longitude
 * (radians) ,height (km) and (Julian) date
 * model=1 is IGRF90, 2 is WMM85, 3 is WMM90, 4 is WMM95,
 * 5 is IGRF95, 6 is WMM2000, 7 is IGRF2000
 * N and E lat and long are positive, S and W negative
*/

  public static double SGMagVar( double lat, double lon, double h, long dat, int model, double[] field ) {
    /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
    int n,m;
    double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
    double sinpsi, cospsi, inv_s;

    boolean been_here = false;

    double sinlat = Math.sin(lat);
    double coslat = Math.cos(lat);
    double mm = 0.0;

    /* convert to geocentric */
    sr = Math.sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
    /* sr is effective radius */
    theta = Math.atan2(coslat * (h*sr + a*a), sinlat * (h*sr + b*b));

    /* theta is geocentric co-latitude */

    r = h*h + 2.0*h * sr +
	    (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
	    (a*a - (a*a - b*b) * sinlat*sinlat );

    r = Math.sqrt(r);

    /* r is geocentric radial distance */
    c = Math.cos(theta);
    s = Math.sin(theta);
      /* protect against zero divide at geographic poles */
    inv_s =  1.0 / (s + ((s == 0.)?1.0e-8:0));

    /*zero out arrays */
    for ( n = 0; n <= nmax; n++ ) {
      for ( m = 0; m <= n; m++ ) {
      	P[n][m] = 0;
	      DP[n][m] = 0;
      }
    }

    /* diagonal elements */
    P[0][0] = 1;
    P[1][1] = s;
    DP[0][0] = 0;
    DP[1][1] = c;
    P[1][0] = c ;
    DP[1][0] = -s;

    /* these values will not change for subsequent function calls */
    if( !been_here ) {
      for ( n = 2; n <= nmax; n++ ) {
      	root[n] = Math.sqrt((2.0*n-1) / (2.0*n));
      }
      for ( m = 0; m <= nmax; m++ ) {
	      mm = m*m;
	      for ( n = Math.max(m + 1, 2); n <= nmax; n++ ) {
	        roots[m][n][0] = Math.sqrt((n-1)*(n-1) - mm);
	        roots[m][n][1] = 1.0 / Math.sqrt( n*n - mm);
	      }
      }
      been_here = true;
    }

    for ( n=2; n <= nmax; n++ ) {
      /*  double root = sqrt((2.0*n-1) / (2.0*n)); */
      P[n][n] = P[n-1][n-1] * s * root[n];
      DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) * root[n];
    }

    /* lower triangle */
    for ( m = 0; m <= nmax; m++ ) {
      /*  double mm = m*m;  */
      for ( n = Math.max(m + 1, 2); n <= nmax; n++ ) {
        /* double root1 = sqrt((n-1)*(n-1) - mm); */
        /* double root2 = 1.0 / sqrt( n*n - mm);  */
        P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
             P[n-2][m] * roots[m][n][0]) * roots[m][n][1];
        DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
              (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) * roots[m][n][1];
      }
    }

    /* compute gnm, hnm at dat */
    switch(model) {
    case 1:	/* IGRF90 */
      yearfrac = (dat - yymmdd_to_julian_days(90,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
          gnm[n][m] = gnm_igrf90[n][m] + yearfrac * gtnm_igrf90[n][m];
          hnm[n][m] = hnm_igrf90[n][m] + yearfrac * htnm_igrf90[n][m];
	      }
      }
      break;

    case 2:	/* WMM85 */
      yearfrac = (dat - yymmdd_to_julian_days(85,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
          gnm[n][m] = gnm_wmm85[n][m] + yearfrac * gtnm_wmm85[n][m];
          hnm[n][m] = hnm_wmm85[n][m] + yearfrac * htnm_wmm85[n][m];
      	}
      }
      break;

    case 3:	/* WMM90 */
      yearfrac = (dat - yymmdd_to_julian_days(90,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
          gnm[n][m] = gnm_wmm90[n][m] + yearfrac * gtnm_wmm90[n][m];
          hnm[n][m] = hnm_wmm90[n][m] + yearfrac * htnm_wmm90[n][m];
        }
      }
      break;

    case 4:	/* WMM95 */
      yearfrac = (dat - yymmdd_to_julian_days(95,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
         gnm[n][m] = gnm_wmm95[n][m] + yearfrac * gtnm_wmm95[n][m];
         hnm[n][m] = hnm_wmm95[n][m] + yearfrac * htnm_wmm95[n][m];
	      }
      }
      break;
    case 5:	/* IGRF95 */
      yearfrac = (dat - yymmdd_to_julian_days(95,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
          gnm[n][m] = gnm_igrf95[n][m] + yearfrac * gtnm_igrf95[n][m];
	        hnm[n][m] = hnm_igrf95[n][m] + yearfrac * htnm_igrf95[n][m];
	      }
      }
      break;
    case 6:      /* WMM2000 */
      yearfrac = (dat - yymmdd_to_julian_days(0,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
           gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
      	   hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
	      }
      }
      break;
    case 7:      /* IGRF2000 */
      yearfrac = (dat - yymmdd_to_julian_days(0,1,1)) / 365.25;
      for (n=1;n<=nmax;n++) {
        for (m = 0;m<=nmax;m++) {
          gnm[n][m] = gnm_igrf2000[n][m] + yearfrac * gtnm_igrf2000[n][m];
          hnm[n][m] = hnm_igrf2000[n][m] + yearfrac * htnm_igrf2000[n][m];
        }
      }
      break;
    }

    /* compute sm (sin(m lon) and cm (cos(m lon)) */
    for (m = 0;m<=nmax;m++) {
      sm[m] = Math.sin(m * lon);
      cm[m] = Math.cos(m * lon);
    }

    /* compute B fields */
    B_r = 0.0;
    B_theta = 0.0;
    B_phi = 0.0;
    fn_0 = r_0/r;
    fn = fn_0 * fn_0;

    double c1_n = 0.0;
    double c2_n = 0.0;
    double c3_n = 0.0;
    double tmp = 0.0;
    for ( n = 1; n <= nmax; n++ ) {
      c1_n=0;
      c2_n=0;
      c3_n=0;
      for ( m = 0; m <= n; m++ ) {
        tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
        c1_n += tmp * P[n][m];
        c2_n += tmp * DP[n][m];
        c3_n +=  m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
      }
      /* fn=pow(r_0/r,n+2.0);   */
      fn *= fn_0;
      B_r += (n + 1) * c1_n * fn;
      B_theta -= c2_n * fn;
      B_phi += c3_n * fn * inv_s;
    }

    /* Find geodetic field components: */
    psi = theta - (pi / 2.0 - lat);
    sinpsi = Math.sin(psi);
    cospsi = Math.cos(psi);
    X = -B_theta * cospsi - B_r * sinpsi;
    Y = B_phi;
    Z = B_theta * sinpsi - B_r * cospsi;

    field[0]=B_r;
    field[1]=B_theta;
    field[2]=B_phi;
    field[3]=X;
    field[4]=Y;
    field[5]=Z;   /* output fields */
    /* find variation in radians */
    /* return zero variation at magnetic pole X=Y=0. */
    /* E is positive */
    return (X != 0.0 || Y != 0.0) ? Math.atan2(Y, X) : (double) 0.0;
  }

  public static double getMag(double lat, double lon) {
//    System.out.println("Starting "+ Magfield.class.getName());

    double height = 0;
    Calendar now = new GregorianCalendar();
    int yy = now.get(Calendar.YEAR) % 100;
    int mm = now.get(Calendar.MONTH) + 1;
    int dd = now.get(Calendar.DAY_OF_MONTH);
    int model = 6;
    double[] field = new double[6];
    int latsign = +1; // N
    int lonsign = +1; // E
    int varsign = +1; // E

//    System.out.println("lat = "+ latsign*Magfield.deg_to_rad(lat));
//    System.out.println("lon = "+ lonsign*Magfield.deg_to_rad(lon));
//    System.out.println("alt = "+ height);
//    System.out.println("julian = "+ Magfield.yymmdd_to_julian_days(yy,mm,dd));
//    System.out.println("model = "+ model);

    double variation = Magfield.rad_to_deg(varsign*Magfield.SGMagVar(latsign*Magfield.deg_to_rad(lat),lonsign*Magfield.deg_to_rad(lon),height,
                Magfield.yymmdd_to_julian_days(yy,mm,dd),model,field));
    
    return variation;
  }
}

